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1) Introduction 

• The lattice Boltzmann equation methodology is a general framework for approximat- 
ing problems that are modelled by partial differential equations under conservative form 
arising in physics. It has been first proposed for fluid dynamics in the context of cellular 
automata (see e.g. Hardy et al [TT], Wolfram |36], d'Humieres et al [6j). The classical 
lattice Boltzmann scheme (McNamara-Zanetti |27) . Higuera et al [20], Qian et al [29], 
d'Humieres has been first developed for nonlinear fluid problems. It has also the 
capability to approximate thermal flows. (Chen a/ [3]), magnetohydrodynamics (Dellar 
[7J) or coupling between fluid and structures (see e.g. Kwon [22j). 

• We have proposed in [Bl E] the Taylor expansion method to analyze formally the 
d'Humieres lattice Boltzmann scheme [5] when the mesh size tends to zero. We have then 
replaced the (much more) formal Chapman- Enskog expansion methodology used in [5] 
by a simple Taylor expansion relative to the grid spacing. In this way we can obtain the 
modified equations of the scheme (in the sense of Lerat-Peyret [2S] and Warming-Hyett 
|34| : see also [TB| [21 |33]) at an arbitrary order in the general nonlinear case. We made 
more explicit in [10] the result at third order accuracy in all generality. We have observed 
at this occasion that a serious difficulty with the lattice Boltzmann scheme lies in the fact 
that the equivalent mass conservation equation contains an a priori non-null third order 
term. We have also proposed an algorithm to derive the modified equation in the case of 
linearized problems [TT]. We have applied this methodology to derive the so-called "quartic 
parameters" to enhance the accuracy of the lattice Boltzmann scheme to simulate shear 
waves [11] . An application of this result is used by Leriche et al [26] for computing with 
spectral method and lattice Boltzmann scheme the eigenmodes of the Stokes problem in a 
cube. This Taylor expansion method can also be used for a detailed analysis of boundary 
conditions. In [12] we have shown that the "magic boundary parameters" of Ginzburg 
and Adler [H] depend in fact upon the detailed choice of the parameters of the lattice 
Boltzmann scheme. 

• In this contribution, we consider linearized athermal acoustics in two and three space 
dimensions. In Section 2, we recall the essential properties of the d'Humieres scheme. 
Then in Section 3, we use the method of formal expansion to expand the discrete eigen- 
modes of the acoustic system of partial differential equations as the mesh size tends 
to zero. In Section 4, we increase the accuracy of the lattice Boltzmann scheme with 
the development of "quartic" parameters and develop also a weaker approach by enforc- 
ing isotropy at fourth order accuracy. Applications in two and three space dimensions 
(D2Q13 and D3Q27 schemes) are presented in Sections 5 and 6 respectively. We detail 
our version of the D3Q27 scheme and the explicit formulae for the determination of the 
quartic parameters. 
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2) d'Humieres lattice Boltzmann scheme 

• In the framework proposed by d'Humieres [5], the lattice Boltzmann scheme uses a 
regular lattice C with typical mesh size Ax and for each node x in £, a discrete set 
of (J+ 1) velocities V is given. In this contribution, the set V does not depend on the 
vertex x. We introduce a velocity scale A and the time step At is obtained according to 
the so-called "acoustic scaling" : 

w - - X 

For X G £ and Vj G V, the point x + VjAt is supposed to be a new vertex of the lattice. 
The unknown of this lattice Boltzmann scheme is the particle distribution fj{x, t) for 
X G £, < j < J and discrete values of time t. The numerical scheme proceeds into two 
major steps. 

• The first step describes the relaxation / — y f* of particle distribution / towards 
a locally defined equilibrium. It is local in space and nonlinear in general. In this paper 
devoted to acoustics we consider only linear contributions. It is defined with the help of a 
fixed invertible matrix M. The moments are defined through a simple linear relation 

J 

(2) mk = J2^kjfj, 0<k<J. 

j=0 

Note that this very simple linear hypothesis is not satisfied in the scheme proposed by 
Geier et al |13| . The first d + 1 moments (where d is the space dimension, equal to d = 2 
or d = 3 in the present applications) are the total density p and the momentum : 



r J 



(3) 



mo = P = ^ 



j=0 

J 



m^ = Qa 



j=0 



We denote by W the vector (in M^) composed of the density and the components of the 
momentum. These moments are supposed to be at equilibrium: 

(4) m* = mi = m";"^ = Wi , 0<i<d. 

The relaxation evolution m — > m* due to linearized collisions is local and trivial for 
the vector W of conserved quantities, as remarked in (jl]). For the other moments, a 
distribution of equilibrium moments G{») is given as a (linear for the present study of 
acoustic waves) function of the vector W: 

(5) ml"" = Gk{W), d<k<J. 
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For k > N, the evolution m — > m* is supposed to be uncoupled: 

(6) ml = mk + Sk {rn^ - rrik) , k > d. 

It is parametrized by the so-called "relaxation rate parameters" Sk- For stability reasons 
(see e.g. Lallemand and Luo we have the inequalities 

(7) < Sk < 2, k> d. 

When the components t) are computed for each x E C at discrete time t, the 

distribution fj{x, t) after relaxation is reconstructed by inversion of relation ([2]): 

J 

(8) /; = 5^M-/m,*, 0<j<J. 

• The second step is just a free advective evolution without collision through charac- 
teristics: 

(9) fj{x,t + At) = f*{x -VjAt,t) , 0<j<J, xeC. 

• The asymptotic analysis of cellular automata provides evidence supporting asymptotic 
partial differential equations with transport coefficients related to the induced parameter 
defined by the so-called Henon's relation |18] 



11 

(10) ak^---. 

The lattice Boltzmann scheme is classically considered as second order accurate (see e.g. 
[23j). We describe in Section 6 the D3Q27 (c? = 3, J = 26) numerical scheme that we 
used for acoustic application. 



3) Formal development of eigenmodes 

• With the method described in [TT], it is possible to derive the equivalent system of 
id + 1) equations associated to a lattice Boltzmann scheme. The system of linearized 
conservation equations issued from this analysis can be written under the form 

(11) A{Ax,d),W = ^ 

which uses a compact notation stating that A(Ax, &) is a (rf + 1) x (c? + 1) (4x4 for 
three-dimensional applications, 3x3 for two-dimensional models) matrix of differential 
operators of high order acting on the conserved variables W . 
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We seek the eigenmodes of the operator A{Ax,d), id est the eigenvalues Xj{Ax,d) 
and the eigenvectors rj{Ax,d) such that 



(12) 



A{Ax, d).rj{Ax, d) = \j{Ax, d) rj{Ax, d) , < j < d. 



We introduce the diagonal matrix A{Ax,d) composed of the eigenvalues Xj{Ax,d) and 
the square matrix R{Ax, d) composed of the eigenvectors. Then relation f lT2|) can be 
written under the following synthetic form: 



(13) 



A{Ax, d) . R{Ax, d) = R{Ax, d) . A(Ax, d) . 



• Moreover, the operator A{Ax, d) in equation fllip is a formal series in the variable 
Ax that is truncated at order 4 in this contribution: 



(14) 



A{Ax, d) = Ao{d) + Ax Ai{d) + Ax^ A2{d) + Ax^ Asid) + 0{Ax^) 



We can apply in this case the so-called stationary perturbation theory for linear operators 
(see e.g. [3j for an elementary introduction). First for Ax = 0, the operator Ao{d) 
is exactly the perfect linear acoustic model. For d = 3 (the case d = 2 is a simple 
consequence of this particular analysis and we omit the details), we have 



(15) 



Ao{d) 



( dt d. dy d,\ 

eld. dt 

eldy dt 

\cld, dtj 



Introducing the Laplacian operator A = d^ + dy + d'^ , a. system Ro{d) of reference 
eigenvectors can be given under the form 



(16) 



Roid) 



/cqVA 


-CqVA 





\ 


Cq dx 


Cq dx 


dy 




Cq Uy 


r^d 

Cq Uy 


-dx 




\cld. 


c^d 

Cq Uz 





dl + dll 



It is elementary to observe the relation 

(17) A,{d).R^{d) = R^{d).Ao{d) 
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and the associated matrix Ao{d) of order zero is simply given by 



(18) 



Ao(5) 



fdt + coVA 0\ 

dt-coVX 

dt 

\ dtj 



We observe that even if the eigenvalue dt is degenerated, i.e. when the eigenvalues are 
multiple, the family of eigenvectors proposed in (ITB]) is complete and defines a basis. This 
allows us to show the existence of two acoustic waves and two shear waves (see e.g. [24J). 

• The parameter Ax is supposed to be infinitesimal and we introduce a formal expan- 
sion of the eigenvalues with diagonal matrices Aj{d): 

(19) A(Ax, (9) =Ao(9) + Aa;Ai(9)+Aa;'A2(9) + Ax='A3(9) + 0(Ax^) 
and relative perturbations Qj{d) of the eigenvectors: 

(20) RiAx, d) = Roid) . (Id + Ax Qiid) + Ax^ Q^id) + Ax^ Qs{d) + 0(Ax^)) . 

We adopt in this work the scaling condition for eigenvectors proposed e.g. by Cohen- 
Tannoudji et al [3J. The component of rj{Ax,d) relative to the corresponding unper- 
turbed eigenvector is exactly equal to unity. In other words, all the diagonal terms of the 
Qj{d) matrices are null: 



(21) 



Q,{d) 



0, 



0<i<d, j>l. 



We insert the expansions fl2Ul) and fl^T]) into the eigenmode condition (IT^ . We introduce 
the expression Aj{d) for the perturbation relative to the unperturbed eigenbasis, id est 

(22) A,id)=Roid)-\A,id).Roid). 

Then we identify each term relative to Ax and we obtain in this manner: 

(23) Ai(9) = Ii(9) + Ao(9) . Qi{d) - Q,{d).Ao{d) 



(24) A2(9) = A2{d) + A,id).Qiid) - gi(9) . Ai(9) + Ao(9) .g2(9) - g2(9) . Ao(9) 



(25) 



A3id) + ^2(9) .Qi(9) - Qi(9) . A2(9) 
A,{d) ^ { +A,{d).Q2{d)-Q2{d).Ar{d) 
+Ao{d).Qsid)-Qs{d).Ao{d) 
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Explicit forms of the eigenvalues at each order may be found without difficulty using 
symbolic manipulation. 

• We introduce a reference density po, a reference sound velocity cq, the shear viscosity 
p and the bulk viscosity (. It is well known (see e.g. Landau-Lifshitz [21]) that the four 
compressible isothermal acoustic linearized equations take the form: 

dtP + Po div u = 
Po dtu + Co d^p - pAu- (C + 5x (div u) = 
podtv + cldyp - pAv - (c + dy{divu) = 



(26) 



Po 



dtw + cld.p- pAw- (^C + ^)M'^'^^) = 



with divu = dxU + dyV + dzW and we recognize matrix ( IT5|) when p = ( = 0. This system 
has two families of eigenvalues. 

• On the one hand, the acoustic waves A± are given by 



(27) A± = - 7A ± coVA Jl + ^A 



2 



c, 



and are parametrized by the attenuation of the sound waves 7 given according to 

We can expand the previous acoustic wave eigenvalues in powers of 7 when the attenuation 
is sufficiently small compared to the frequency of the sound waves. Then we have the 
expansion 

(29) X± = dt ± co^/A-^A ±l^A'/^ + Oi^""). 

2 Co 

The previous (relatively unusual) notations with matrices of operators are exactly the 
one used when implementing without difficulty (due to linearity) our approach with 
a symbolic manipulation software. As is well known, the pseudo-differential operators 
like a/a and A'^/^ are defined by their action on Fourier transforms (we refer e.g. to 
the book of Hormander |2T]). It is also possible to introduce a Fourier decomposition 
on harmonic waves of the type exp(i(a;t — k.x)). Then we have the usual change of 
notation: dt = ico, V = — ik, A = — |kp, A'^/^ = —i |k|^, etc. The corresponding 



dispersion relation associated with the annulation of the eigenvalues \± given in ([29 
allows the introduction of the parameter according to 

2 

(30) -ico = Te = J |k|2 ± tco |k| (l-^ IM^) + 0(7^ |k|^) . 
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• On the other hand, the shear wave Aq define a multiple eigenmode that is simply 
given by 

(31) Xo = dt-jyA 

with u = —. With the spectral notation, relation (1311) can be written with the help of 
the attenuation Fj of the shear waves under the form 

(32) -iu = Tt = 1^ |k|2 . 

As previously, a family of two independent eigenvectors with two different "polarizations" 
corresponds to this shear eigenvalue. Moreover, the perturbation Ax Ai{d) + Ax'^ ^2(8) + 
Ax^ A^ld) + ■ ■ ■ of relation f ll4p is diagonal on the basis of unperturbed eigenvectors 
proposed at relation (fT6|) . Due to well-established arguments (see e.g. the textbook of 
Quantum Mechanics of Cohen- Tannoudji et al [3j or e.g. the sections "multiple roots" and 
"degenerate roots" in the book of Hinch [H]), the expansion (120|) f l2T]) is completely valid 
even in this (relatively) complicated degenerate case. 



4) Increasing the accuracy of shear and acoustic waves 

• When we use a lattice Boltzmann scheme, the viscosity admits a representation of 
the type (c./. the D3Q27 scheme): 

(33) /i = - A Ax 

and similarly the bulk viscosity is given for the D3Q27 scheme (to fix the ideas) ac- 
cording to 

(34) C = AAxae(^-c^ 

with the notations of Section 6. We can identify the developments (fT9|) on the one hand 
and (l29|) f l3T]) on the other hand because, due to (133|) . f lM|) and (128|) . the viscosities 7 
and IX are proportional to Ax. The usual implementation of lattice Boltzmann scheme 
consists in adjusting the eigenvalues Ao(i9) and Ai((9) in order to fit the viscosities. 

• When we identify the developments (1191) and 0311) at the order 4 for the shear waves, 
we recover the "quartic parameters" obtained previously [Hj, in particular for the D2Q9 
[23] and D3Q19 [32j schemes. In these cases, the shear viscosity ^ follows the "quartic" 
condition 

(35) u = ^^AAx. 
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If we add to the previous study the fitting of the acoustic waves ( 129|) . there is no solution 
of the equations for the two previous schemes D2Q9 and D3Q19. We have then considered 
the previous conditions for extended stencils such as the D2Q13 PU| 155] and D3Q27 (see 
e.g. [28j) lattice Boltzmann schemes. The results for quartic coefficients are displayed in 
Section 6 when fitting all the waves (129|) and f l3T]) of the linear problem. 

• In order to take advantage of the fiexibility of the d'Humieres scheme, we have devel- 
oped an "isotropic" condition, not as strong as the quartic condition. With this isotropic 
condition, we do not impose anymore the conditions that the numerical waves (IT^ fit ex- 
actly the physical ones (129|) and f l3T]) but suggest that the numerical waves are isotropic. 
In other words, the operator matrices A2{d) (for the order 3) and A3 (9) (for the order 4) 
are differential operators that are constrained in order to depend only on the radial vari- 
able r (r^ = x"^ +y'^ + z'^) and on the operator d/dr . When this condition is imposed, it is 
possible to fit isotropic waves for the four previous schemes. In Table 1 and 2, we display 
the number of (a priori non-independent) discrete equations that must be satisfied by 
the parameters of the lattice Boltzmann schemes when we impose isotropy for the shear 
waves, isotropy for the acoustic waves, exact fitting of the shear waves and exact fitting 
of the acoustic waves. All these discrete equations have been solved exactly with the help 
of symbolic manipulation. To fix the ideas, we give in Section 6 the parametrization of 
the quartic shear and acoustic waves for the D3Q27 lattice Boltzmann scheme (solution 
of 18 discrete equations, according to Table 2). 





isotropic shear 


isotropic acoustic 


exact shear 


exact acoustic 


order 3 





1 





2 


order 4 


1 


1 


2 


2 



Table 1. Fitting "isotropic" and "exact" shear and acoustic waves at various orders of 
accuracy. Number of equations for the D2Q9 and D2Q13 lattice Boltzmann schemes. 
Note that the same numbers apply to both lattices. 





isotropic shear 


isotropic acoustic 


exact shear 


exact acoustic 


order 3 





1 





2 


order 4 


10 


2 


13 


3 



Table 2. Fitting "isotropic" and "exact" shear and acoustic waves at various orders of 
accuracy. Number of equations for the D3Q19 and D3Q27 lattice Boltzmann schemes. 
Note that the same numbers apply to both lattices. 
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5) Linearized acoustics with the D2Q13 scheme 

• Two numerical "experiments" have been performed: emission from the center (Figures 
1, 2 and 3) and emission from the boundary of the circle (see Figures 1, 4, 5, 6). In both 
cases, a first order anti-bounce-back boundary condition (see e.g. [IJ) is implemented 
to impose a given density (pressure) on the boundary: either a constant when sound 
is emitted at the center, or a sinusoidal time dependent value for the emission by the 
boundary. In each case, we measure the value of the density at a fixed position in space 
vs time (see Figure 6) or we analyze the pressure field at a given time (see Figures 2 
to 5). Obviously the maximum duration of the simulation with a radius R and sound 
velocity cq is R/cq when the source is at the center and 2R/co when it is on the circular 
boundary. These experiments are sensitive to anisotropic behaviour of both the velocity 
and of the attenuation as can be seen in the figures. 





Figure 2. D2Q13 Acoustic propagation, diverging acoustic test case, time period = 8 
(7 points per wave length, Cq = O.SOj.' usual (left), isotropic (center) and quartic 
(right). 
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-0.015 ' ' ' ' ' 1 -0.02 ' ' ' ' ' 1 

85 90 95 100 105 110 85 90 95 100 105 110 

radius radius 



Figure 3. D2Q13 diverging acoustic test case. Measurements after 200 time steps; 
usual parameters (left) and quartic parameters (right). 




Figure 4. D2Q13 Converging acoustic propagation, time period = 8 (7 points per wave 
length, Cg = O.SOj; usual (left), isotropic (center) and quartic (right). 




Figure 5. D2Q13 Converging acoustic propagation, time period = 10 (9 points per 
wave length, Cq = 0.80j; usual (left), isotropic (center) and quartic (right). 
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240 260 280 300 320 340 360 330 400 240 260 280 300 320 340 360 380 400 

time time 



Figure 6. D2Q13 Converging acoustic propagation, time period = 8 (left) and time 
period = 10 (right). 



6) D3Q27 athermal linearized Acoustics 

• The D3Q27 Lattice Boltzmann sdieme, described with details e.g. in Mei et al. [2S] 
is simply obtained by taking a tensorial product for the set V of discrete velocities: 



(36) 



^3 



-1,0,+1, 0<j<26. 



Due to the large number of moments, we detail in this subsection the way the matrix M 
parametrizing the transformation (j2]) is obtained. First, velocities f " for < j < J = 26 
and 1 < a < 3 are given according to relation ( 136|) . The first four conserved moments 
p and are determined according to ([3]). The generation of other moments is done 
according to the tensorial nature of the variety of moments that can be constructed, as 
analyzed by Rubinstein and Luo jSI]: scalar fields are naturally coupled with one another, 
idem for vector fields, and so on. So components of kinetic energy are introduced: 



(37) 



m4 



26 3 



j=0 13=1 



The entire set of second order tensors is completed according to 



(3J 



26 



j=0 
26 



j=0 
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(39) 



26 



j=0 
26 

26 
j=0 



Flux of the energy and of the square of energy: 

26 3 



(40) 



Square and cube of the energy: 



3 E(E 

j=0 /9=1 
„ 26 3 



1< a < 3 



(41) 



26 3 



mi6 = e 



mn = eg = 



j=0 /3=1 
„ 26 3 „ 

E(E(''f)1 



3=0 /3=1 



Product of XX and WW by the energy: 



(42) 



26 



m^s -XX^^ 3 (2 iv]r - {v]r - ^ (.f)^) /. 

3=0 /3=1 
26 3 

m„ = H-'W. ^ 3 ^ ((.J)^ - (vff) (E ("i")') * 



J=0 



Extra-diagonal second order moments of the energy: two different components times 
energy and permutations 



(43) 



3=0 13=1 

26 3 

3=0 /3=1 

26 3 



i=o 



/3=1 



14 



Frangois Dubois and Pierre Lallemand 



Third order pseudovector 

26 

(44) m22+. = = E < - (^r')') A ' 1 < « < 3 . 

i=o 

with a natural permutation convention and the third order totally antisymmetric tensor 
XYZ: 

26 

(45) m26 = XYZ = ^ v] v] f, . 

j=0 

A first (non-orthogonal) matrix M is obtained by using ([2]) ([27]) ([M]) (El]) (110]) (111]) (i^) 
( H3|) dH]) and (145|) as linear combinations of the //s : 

26 

(46) mk = J2 Mkj fj , 0<k<26. 

j=0 

Then matrix M is orthogonalized from relations (jS]), (EZD, (jSH]), §9^, (gO]), (gl]), (H2|). 
f H3|) . (jH]) and (H5|) with a Gram-Schmidt orthogonalization algorithm: 

Mij = Mij - 9ie Mej, z > 4 . 
The coefficients ^^j^ are computed recursively in order to enforce orthogonality: 

26 

Mij Mkj = for i^k. 

j=0 

Note that the coefficients in ([37D, ([38D, (gO]), dH]), (12]), (gS]), (gH) and (g5]) have 

been chosen such that the coefficients of matrix M are all integers. 

• The equilibrium values of the moments define the equilibrium distribution. The first 
moments p and are in equilibrium and have respectively a scalar and a vectorial 
structure. The equilibrium values of the other moments respect this structure. We have 





= ex'p 












= WW^'i = 




= = 


ZX^^ 


= 




= Ci Qa 
























= /3AV 










^3 


= eAV 












= WW^'i = 




= = 


ZX|<^ 


= 


a 


= C3 A^ Qa 










XYZ'''' 


= 0. 











V 
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We have the following relation between the parameter 6 and the sound velocity cq 



(4J 



9 = 3cl-2 



The relaxation parameters for these moments are denoted respectively by the following 
values 



(49) 



xx'''^, ww'"^, zx""^ 

s-eq 



xxi"^, WW^ 

XY^, YZl"^, ZX^ 

a 

XYZ^'^ 



Se 

St 
Suj . 



With these relaxation parameters, the shear viscosity is isotropic when 



(50) 



ci = -2 

C3 = 



• The quartic parameters for athermal acoustics are solution of the 13 + 2 + 3= 18 
discrete equations obtained at Table 2. They are displayed in this subsection. 



(51) 



1 

2 {a^ - cxe) 



9 cl - 18 cl ae + 27 at + 180 cla^ a, 
+144 clal-8a^ + 8ae- 324 a^c (^e) 



(52) 



C2 = -42(7^ + - 
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(54) <^ 



_1 ^ 
48 D, 

-76 al + 27 cl- 27 4 

-93312 4 at a'^. 



180 4 al - 468 4 (^1 + 324 4 al + 7776 4 (^1 



4 at + 80 a^ a^ - 336 at a^ 



13AAalal + 240at 



10800 4 al (Te - 46656 4 (^1 (^l + 20736 al al 4 + 62208 4 al a, 



3 „3 /I 



-324 de + 3888 cx^ al + 864 al - 4752 al 
+51840 d^. + 3888 4 ^l - 46656 4 ^\ + 324 4 a^ a^ 



-3456 al + 2880 al + 20736 eg + 8064 al al + 6912 al a^ 



-864 eg al + 1440 



14400 al al + 324 4 al + 31104 eg al al 
„ .f^x - f^e) ( - 84 al + 432 eg al + 84 + 540 eg a^. - 23 a^ 
-972 e^ a^ al - 27 eg a,. + 81 e^ a,. ^23 a,- 54 eg ) 



a. 



We observe that the presence of the factor [ox — o"e) in the denominator of the algebraic 
expression of a^ makes this expression incompatible with a BGK or TRT [ISj type 
hypothesis. The same remark holds for relations (155|) and ( 156|) relative to the parameters 
(T^ and cr^ respectively. 



(55) 



a^ 



1 



336a,(-l + 12a2.) (a,-(Te)2 



1968 a^- 144 eg cr^ 



15552 4 al 



mAalal 



1344 < ae 



+4608 ae + 8064 al al 

A 4 2 



12672 a° + 103680 e^a>, 

186624 e^cT>" 
+80 (7^ ae -76al + 82944 eg + 15552 4 al 



27 4 + 274-4: al + 324 4 a] - 180 eg al 



^ 96a,(-l + 84a2)(cr,-ae)2 ^ 

A^;^ = 192 al - 828 eg al + 972 e^ + 6480 eg al al - 2592 eg al 

-11664 4 al al + 192 al al - 384 cr^ a^ + 2304 cXe + 4032 al al 
-6336 al + 7Qal + 51840 eg al - 93312 e^ al + 27 e^ - 27 eg 
+4 al - 324 e^ + 180 eg al - 80 a,, a^ + 41472 eg al + 7776 e^ al 



(57) 



1 Nr 

12 'Or 

76a^ - 324e^ae^ + 180ega^ 
+720 al al 



144 al - 80 a^ (Je + 3456 4 a 



2 _4 



576 al ae - 504 eg al + 4320 eg al < + 27 e^ + 4 cr^' 



-27eg- 7776e^a,V^ + 648e^a^ 
Dr = ax {7Q al - 528 al - 50A 4 al + QA8 4 al + A320 4 alal 

+3456 eg al - 7776 e^ al + 336 a,^ (x^ + 192 al ae + 27 e^ - 27 eg 
+4 cTe^ - 324 4 (^l + 180 eg al - 80 (x^ ae) 
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(5J 



1 ^ 
12 A, 



816 - 828 eg cr^ + 972 + 6480 a,^ - 2592 a. 



2 _2 2 



„2 _4 



-11664 ci^crfa,^ 



+3456 gI + lQat + 51840 al a'i - 93312 ol a't + 21cl- 27 Cq 
+4 cr,2 - 324 a] + 180 - 80 + 41472 + 7776 
a, ( - 192 al - 828 eg al + 972 + 6480 eg cr^ cr,^ - 2592 eg 
-11664 ct al al - 192 al al + 384 ae + 2304 ae + 4032 al al 



ou^u^- ul6 < < + 1632 < de - 17280 < ae + 13824 a, 

2 I r-io^n„2 4_2 riooio„4 4 2 



4 ^2 

,2 



-6336 al + 76 + 51840 eg a^ al - 93312 e^ a^ al + 27 e^ - 27 eg 
+4 al - 324 e^ al + 180 eg - 80 a^ a^ + 41472 eg a^ + 7776 e^ ) 

For the numerical experiments described above, we have used the following values of the 
equilibrium parameters. We did the previous theoretical numerical computations f lST]) 



software: 



(59) 



Se 
Co 



55|) ( 156|) f l57j) fl58|) with a 50 Digits option in our symbolic manipulation 

= 0.95057034220532319391634980988593155893536121673004 

= 1.8552875695732838589981447124304267161410018552876 

= 0.623538 

= 2.436118000 

= 1 



and the "quartic" relaxation parameters: 



(60) 



/3 

Se 
St 



0.50345521670787922851706691010021052631578947368420 
0.37925445705024311183144246353322528363047001620745 
1.3 

0.34253657030513141274711235609982461596733955718034 
1.2 

1.9945477114942149093456286590460711496091258797386 
1.2940799466197218037166471960307116873345869400515 
1.9451616927239606019667013153498030552793202566039 
0.25131560984615404581501005329813711811837082814760, 



Of course, we used the previous numbers with the usual truncation allowed by the 32 bit 
or 64 bit floating point arithmetic. We observe that for these particular parameters, due 
to the link (155]) between the shear viscosity and the parameter a^-, to relation f lM|) 
between the bulk viscosity and the parameter ae and to formula (128|) making explicit 
the attenuation 7 of sound waves, we have 



(61) 



0.013 A Ax 
0.0920492 A Ax 
0.0546913 A Ax, 
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• The results of simple periodic waves are displayed in Figures 7, 8 and 9. The numerical 
results confirm that parameters proposed in (l59i) and (160|) allow one to get fourth order 
(relative) accuracy. 
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Figure 7. Periodic wave with the lattice Boltzmann scheme D3Q27. Error for shear 
eigenvalue Tt defined in (E^- For quartic parameters, we have — z/ |k|^= 0(|k|^). 
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Figure 8. Periodic wave with the lattice Boltzmann scheme D3Q27. Error for the 
imaginary part of the acoustic eigenvalue Ti defined in / T^Oj) . For quartic parameters, 
we have fifth order accuracy: \ni(Vi^ — Cq |k| (l — 7^ |kp /(2cq)) = 0(|k|^) with 7 
detailed at relation (f^j. 
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Figure 9. Periodic wave with the lattice Boltzmann scheme D3Q27. Error for the real 
part of acoustic eigenvalue Ti. For quartic parameters, we have Ile(Tej — 7 |k|^ = 

o(ikr). 



The numerical experiments have been done for a three-dimensional converging acoustic 
wave. A pulsating sphere of radius i? = 46.08 is embedded in a 95^ cube. The simulations 
used only 64 bit or even 32 bit data with CUDA: (using a Nvidia 9800GT card, one full 
update of D3Q27 takes approximately 11 nanoseconds per site). The "linear anti-bounce- 
back" numerical boundary conditions (in the sense given in Bouzidi et al. [Ij) are used 
to impose a boundary density on the sphere oscillating with a period T = 10. Since 
the sound velocity is around 0.62 (see (15^ ). there are around six mesh points by wave 
length. In Figure 10, the scatter plot of density vs radius after 82 iterations is displayed. 
The results are of excellent quality with the use of quartic parameters. Nevertheless, 
since the attenuation of the sound wave is relatively important for the above parameters 
(7 ^ 0.0547 as recalled in (EI])), the range of propagation is relatively limited (five wave 
lengths in our case, see Figure 10). For the less stringent "isotropic" case, i.e. to fix the 
ideas 



(62) 



( 


1 




6 (Tj. 


< /3 


= 4-9c2 




1 


Co 


7! 



we can find situations of "TRT type" satisfying 

I o"e = (T^ = = = cr^ = 



(63) 
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with much lower attenuation of sound. An example is given in figure 10 with parameters 
of the scheme following relations (162|) and ( 163|) with cr^ = 0.006, that is shear viscosity 
0.002 and sound attenuation 0.002. These results show that choosing the goal in terms 
of accuracy or small attenuation will infiuence the choice of the parameters in the lattice 
Boltzmann simulation. 




5 10 15 20 25 30 35 40 45 

distance from the center 



Figure 10. Sound wave emitted by a sphere. Numerical results with the D3Q27 lattice 
Boltzmann scheme with "TRT" isotropic parameters described at relations / f^) 
and with quartic parameters given by relations / f5P|) and / f^) . Six points by wave 
length. 



7) Conclusion 

• We have proposed to use the Taylor expansion method in conjunction with the use 
of symbolic manipulation to increase the accuracy of the lattice Boltzmann scheme in the 
multiple time relaxation approach proposed by D. d'Humieres for linear acoustic waves. 
The result is the use of the previous scheme with very particular "quartic" parameters. We 
have obtained a family of such parameters for the D3Q27 numerical scheme. Numerical 
experiments confirm the predictions of the theoretical analysis. The catch is to make sure 
that situations are stable. This problem cannot be solved in practice with developments 
around small wave numbers. In the absence of efficient techniques to find stable situations 
with small viscosities and/or sound attenuation, we tried at random a number of sets of 
the free parameters and have used the best of them for the explicit results shown above. 
We note that the less stringent "isotropic" condition is compatible with small viscosities 
which may be quite useful for practical simulations. Comparison of various stencils is 
under way and will be presented elsewhere. 
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